Effective speed of sound in phononic crystals 
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A new formula for the effective quasistatic speed of sound c in 2D and 3D periodic materials is 
reported. The approach uses a monodromy-matrix operator to enable direct integration in one of 
the coordinates and exponentially fast convergence in others. As a result, the solution for c has a 
more closed form than previous formulas. It significantly improves the efficiency and accuracy of 
evaluating c for high-contrast composites as demonstrated by a 2D example with extreme behavior. 
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I. INTRODUCTION 

Long-standing interest in modelling effective elastic 
properties of composites with microstructure has sub- 
stantially intensified with the emerging possibility of de- 
signing periodic structures in air^ and in solids^ to form 
phononic crystals and other exotic metamaterials, which 
open up exciting application prospects ranging from neg- 
ative index lenses to small scale multiband phononic 
devices^. This new prospective brings about the need for 
fast and accurate computational schemes to test ideas in 
silico. The most common numerical tool is the Fourier 
or plane- wave expansion method (PWE). It is widely 
used for calculating various spectral parameters includ- 
ing the effective quasistatic speed of sound in acoustic 4 
and elastic^ phononic crystals. At the same time, the 
PWE calculation is known to face problems when ap- 
plied to high-contrast composites^, which are of especial 
interest for applications. Particularly riveting is the case 
where a soft ingredient is embedded in a way breaking 
the connectivity of densely packed regions of stiff ingre- 
dient. Physically speaking, the speed of sound, which is 
large in a homogeneously stiff medium, should fall dra- 
matically when even a small amount of soft component 
forms a 'quasi-insulating network'. Note that this case, 
which implies a strong effect of multiple interact ions , is 
also ungainly for the multiple-scattering approach^!. 

The purpose of present Letter is to highlight a new 
method for evaluating the quasistatic effective sound 
speed c in 2D and 3D phononic crystals. The idea is 
to recast the wave equation as a lst-order 'ordinary' dif- 
ferential system (ODS) with respect to one coordinate 
(say X\) and to use a monodromy-matrix operator de- 
fined as a multiplicative (or path) integral in x\. By 
this means, we derive a formula for c whose essential 
advantages are an explicit integration in x\ and an expo- 
nentially small error of truncation in other coordinate(s). 
Both these features of the analytical result are shown to 
significantly improve the efficiency and accuracy of its 
numerical implementation in comparison with the con- 
ventional PWE calculation, which is demonstrated for 
a 2D steel/epoxy square lattice. The power of the new 
approach is especially apparent at high concentration / 



of steel inclusions, where the effective speed c displays a 
steep, near vertical, dependence for / w 1, a feature not 
captured by conventional techniques like PWE. 



II. EFFECTIVE SPEED: 2D ACOUSTIC WAVES 

A. SETUP. Consider the scalar wave equation 

V • (fiVv) = -pu 2 v, (1) 



for time-harmonic shear displacement w(x, t) = w(x)e _IW * 
in a 2D solid continuunP with T-periodic density p(x) 
and shear coefficient /i(x). Assume a square unit cell T = 
^2iti&i = [0,1] with unit translation vectors ai _L a 2 
taken as the basis for x = J^. a;^. Imposing the Floquet 
condition v(x) = w(x)e lk ' x where m(x) is periodic and 
k = fc/c ( | ac | = 1), Eq. ([!]) becomes 

(Co + Ci + C2)u = pu 2 u with Cqu = — V(/iVu), 
C\u — — ik • (/iVu + V(/iu)), C2U = k 2 fiu. (2) 

Regular perturbation theory applied to |2]) yields the ef- 
fective speed c(k) = lim w .fc_>o ui(k)/k in the forrrP 

c 2 (k) = /U e ff(K)/(p), Mcff(rc) = (m) - M i K ) witri ( 3 ) 

M ( K ) = I2lj=l M ij K i K 3' M ij = ( C Q ld i^9jfj) = Mji, 

where di = d/dxi, spatial averages are defined by 

(/)=/ T /(x)dx (=((/> 1 ) 2 , (f)i = ftfWdxi), (4) 

and (•,•) denotes the scalar product in L 2 (T) so that 
(/, h) = (fh*) (* means complex conjugation). The dif- 
ficulty with ([3| is that it involves the inverse of a par- 
tial differential operator C - O ne solution is to apply a 
double Fourier expansion to Cq 1 and <9;/i in |3j). This 
leads to the PWE formula for the effective speech which 
is expressed via infinite vectors and the inverse of the 
infinite matrix of Fourier coefficients of /i(x). Numeri- 
cal implementation of the PWE formula requires dealing 
with large dense matrices, especially in the case of high- 
contrast c omp osites for which the PWE convergence is 
slow (see [IV I. An alternative "brute force" procedure of 
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the scaling approach is to numerically solve the partial 
differential equation Coh = dip for the 1-periodic func- 
tion /i(x) (e.g. via the boundary integral methocP). 

The new approach proposed here leads to a more effi- 
cient formula for c based on direct analytical integration 
in one coordinate direction. There are two ways of doing 
so. The first proceeds from the ODS form of the wave 
equation ([lj itself, which means 'skipping' This is 
convenient for deriving c(/c) in the principal directions 
n || ai, 2 , see §IIB. The second method is more closely 
related to the conventional PWE and scaling approaches 
in that it also starts from ([3]) but treats it differently, 
namely, the equation C$h — a-// is cast in ODS form and 
analytically integrated in one coordinate. This is basi- 
cally equivalent to the former method, but enables an 
easier derivation of the off-diagonal component M12 for 
the anisotropic case, see §IIC. 

B. Wave speed in the principal directions. The 

wave equation may be recast as 



= Qv 


A- puj 2 



with A 

,-1 







= - d 2 {pd 2 ) 
jj(x) 



V 

p.v' 



(5) 



where / stands for d\. The solution to Eq. ^ for initial 
data jj(0, x 2 ) = 17(0, •) at x\ = is 



ri(xi, •) = M [xi,0]jj(0, •) with 
M[a,b] = ] a b (I+Qdx 1 ). 



(6) 



The operator M [x±, 0] is formally the matricant, or prop- 
agator, of ([5]) defined through the multiplicative integral 

J (with I denoting the identity operator) . It is assumed 
for the moment that p(x) and /i(x) are smooth to en- 
sure the existence of M. The matricant over a period, 
M [1,0], is called the monodromy matrix. 

Assume the Floquet condition with the wave vector 
k = (ki 0) T so that v(x) = u(x)e ifclXl and 17(1, •) = 
77(0, ■)e lkl . By Mi, this implies the eigenproblem 



M [l,0]w(fei) = e lfcl w(fc 1 ) 



(7) 



where M [1, 0] depends on to. Eq. ^ defines fci = ki(uA 
and hence to — u(ki), where to 2 is the eigenvalue of M 
with i?(x) = tt(x)e The effective speed c(ki) = 

linia>,fci->-o ui/ki can therefore be determined by applying 
perturbation theory to ([7| as ui, ki — > 0. The asymptotic 
form of M [1,0] follows from definitions ([5| and ^ 2 as 

M [1, 0] = Mo + uj 2 Mx + 0(uj 4 ) where : 
Mo = Mo [1, 0] , Mo [a, b] = J°(X + Q d Xl ) with 

'0 /i- 1 



Qo = Qu 



A 

N 



(8) 



Mi = J .Mo [1,2:1] Mo k.i..O] <!./■,. 



Note the identities Qow = 0, Qq w = and hence 

M [a, b] w = w , Mo [a, b] w = w (Va, b) 

for w = (1 0) T , w = (0 1) T . (9) 

By j9ji w is an eigenvector of Mo with the eigenvalue 1, 
and it can be shown to be a single eigenvector. Therefore 
w(fei) = w + fciWi+A;fw2-l-O(fci) and u = cki + 0(k 2 ). 
Insert these expansions along with Qi in ([7]) and collect 
the first-order terms in ki to obtain 



Mo^i = Wi + «wq 



Wi 



{{Mo-1)- 1 ^. (10) 



According to ([9]), Mo — I has no inverse but is a one-to- 
one mapping from the subspace orthogonal to Wo onto 
the subspace orthogonal Wo; hence, Wi exists and wq-Wi 
is uniquely defined. The terms of second-order in ki in 
([7]) then imply 



M0W2 + c .M1W0 = iwi + W2. 



(11) 



Scalar multiplication on both sides by Wp leads, with ac- 



count for 



by (10,) 



Litrpiication 
(|[ and g 



4, to c 2 (p) = —i (wq • Wi) 2 , whence 



c 2 ( Ki )= (p)" 1 (w -(Mo -^-Vo)^ (12) 

where the notation (-) 2 is explained in Q. Interchanging 
variables xi x 2 in the above derivation yields a similar 
result for c(k 2 ) as follows 



c 2 (k 2 ) = (p) 1 (w • (M - I) Vq) where 



So 



Mo = J Q {X + Qodx 2 ) 
'0 pr v 

yA 



(13) 



A 



dindi 



[0,T\ 



The result for a rectangular lattice with T 
[0, T2] is obtained by replacing Xi with Xi/Ti. 

C. The full matrix Mij. The anisotropy of the 
effective speed c(/c), i.e. its dependence on the wave 
normal k = k/fc, is determined by_the quadratic form 



M(k) — J2i j=i MijKiKj (see Eq. «3h), and represented 



by the ellipse of (squared) slowness c (k). Eqs. (12) 



and |13ji, which define c(ni) and so Ma, suffice for 
the case where T is rectangular and /i(x) is even in (at 
least) one of Xi so that the effective-slowness ellipse is 



'(Ki)K~ with the principal axes par- 



allel to ai _L a2. Otherwise c(k) for arbitrary k requires 
finding the off-diagonal component M12. For this pur- 
pose, with reference to ([3]), consider the equation 



Coh = dip: 



(14) 



for 1-periodic h(x). With the above notations this can 
be written as — (fih')' + Ah = p! or, more conveniently, 
[ph'Y = Ah with h — h + xi. The latter is equivalent to 



£' = Qo£ where £ = 



h + xi 
p(h' + 1) 



(15) 
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and Qo is given in (18]) 3. The general solution to |l5|) is 



M o [ Xl ,0]t(0,-), 



(16) 



where Mo [xi,0] is defined in ^2, and£(0, •) is the initial 
data at X\ = 0. The periodicity of h implies £(1,-) — 
£((},•) + w , while = MaW,-) by Hence 

£(0, •) = (Mo -l)~ 1 w and so (14 1 is solved by 



Z( Xl ,-) =M [x 1 ,Q}(Mo-l)- 1 w - (17) 
Substituting ( p~7| ) into the definition of Mi 2 in ^ yields 

M 12 = (C ~ 1 difj,,d 2 fx) = (hd 2 fj) = (9 2 /iw •£) 

= (a 2M w • Af Ni,0](M -2) _1 w ) ■ (18) 



Note that the formula ( 18 ) for M12 requires more com- 
putation than the formulas ( 12 ) and ( 13 )i for Ma. Inter- 



estingly, if the unit cell T is square, then, for an arbitrary 
(periodic) i«(x), Eq. ( |18| ) can be circumvented by using 
the identity M 12 = (M u — M 22 )/2, where My follow from 
Eqs. ( 12 1 and ( 13 )i applied to the square lattice obtained 



from the given one by turning it 45°. 

D. Discussion. The two lines of attack outlinedmen- 
tioned in §11. A are equivalent in that the formula ( 12 ) 



for the effective speed c(k\) in the principal direction 
can also be inferred from Eq. ([3|. Inserting the solution 
(171 of (fl4|) defines the component M\\ as 



M n = (Co^M = (hp!) = (//wq-0 -<zi//)- (19) 

Integrating by parts each term in the last identity and 
using the periodicity of //(x) along with Eqs. pL, ([9]), 
([T5|)-(fT7|) (see also the notation Q) yields 



(Ao ■ = -( w o • (Mo - 1) l vfa)2 + (m(0, X 2 ))2, 
- (xin') = (/i> - (Ai(l,ar 2 )) a - (ji) - (v(0,x 2 )) 2 . (20) 

Thus, Mn= (/i)— (wq ■ (Mo — Z)~ 1 w ) 2 which leads to 



(12), QED. Note that Eq. (18) is also obtainable via the 
monodromy matrix of the wave equation (fTl) (the ap- 
proach of §IIB) with v(x) = it(x)e lk ' x and k jfa^, but this 
method of derivation of M12 is lengthier than in §IIC. 

As another remark, it is instructive to recover a known 
result for the case where //(x) is periodic in one coordi- 
nate and does not depend on the other, say p(xi,x 2 ) = 
H(xt). Using d8b 2 , ®3 and JT3} 3 gives 



(M -l) 





1 



= w , (Mo -I) 





/i(xi 



= w . 



Therefore, by ( 12 ) and (131 1, c 2 («i) 



(21) 



" *>iV<p) and 



c z (n 2 ) = Wi/M while M X2 = by Q with <9 2 /^ = 0. 

Finally, we note that, while the above evaluation of 
quasistatic speed c is exact, using the same monodromy- 
matrix approach also provides a closed-form approxima- 
tion of c. For the isotropic case, it is as follows (see^ for 
more details): 



1 



2(P) 



(^ 1 )7 1 ) 2 + (^) 5 



(22) 



III. EFFECTIVE SPEEDS IN PRINCIPAL 
DIRECTIONS FOR 3D ELASTIC WAVES 

The equation for time-harmonic elastic wave motion 
v(x, t)=v(x.)e~ zu ' t is, with repeated suffices summed, 

- dj(ci jkl div k ) = poj 2 Vi (i,j, k,l = l, 2, 3), (23) 

where density p(x) and compliances cyfci(x) are T- 
periodic in a 3D periodic medium. Assume a cubic unit 
cell T = ti a i — [0> l] 3 and refer the components Xi, v.i 
and Cijki to the orthogonal basis formed by the transla- 
tion vectors a^. Impose the condition v(x) = u(x)e lk x 
with periodic u(x) = (ui) and take k parallel to one of 



aj, e.g. to a!. Eq. (23) may be rewritten in the form 
V ' = Q V with »j(x) = (, ^ ^ , 

' \{CilklOlUk) J 
e=f..2„X rT^A A A C 'Zi) (24) 



v w 2 p<% +A 2 - AxC^Ai AxC- 1 ) 

where the self-adjoint matrix operators C and A\ t2 are 

C = (cfifci)) Ai(ui) = (c ilka d a u k ), 
A 2 (ui) = (d a (c iakb d b u k )) with a, b = 2, 3. (25) 

Like in the 2D case, denote the monodromy matrix for 

(24) at u = by Mo = J Q (l+Q Q dx x ) where Q = Q w=0 , 
and also introduce the 6x3 matrices Wo = (<% 0) T and 

Wo = (0 Sij) T . Reasoning similar to that in §11. C leads 
us to the conclusion that the effective speeds c a (n\) = 
linv/d^o w/fci (a — 1, 2, 3) of the three waves with k = 
kn parallel to ai are the eigenvalues of the 3x3 matrix 

= &)■ (26) 



W • (M -X)- 1 W )J 3 (willi 



IV. NUMERICAL IMPLEMENTATION 

There are several ways to use the above analytical re- 
sults for calculating the effective speed. One approach 
is to transform to Fourier space with respect to coordi- 
nate^) other than the coordinate of integration in the 
monodromy matrix. Consider the 2D case and apply the 



in x 2 



Fourier expansion f(xi,x 2 ) = Y^ n & fn(^i)e 
for the functions / = // and [iT 1 . Then the operator of 
multiplying by the function p~ 1 (xi, •) and the differential 
operator A(x\)= —d 2 (p(xi, -)d 2 ) become matrices 



(Pn—rn) , 



/j 1 1 — > fi 1 (x 1 ) = (p, 1 n _„ 
A 1 — > A(xi) = Air 2 (nrrij2 n - m ) , n^m^TL, (27) 



and Eq. ( 12 ) reduces to following form 
c 2 (K 1 ) = {p)- 1 w d -(M 

M = J^I + Qodxi) 

= (0 S 0n ) T , wg = (Son 0) T 



I) 1 Wq with 

Qo(xO = (l ^ 



(28) 
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where c{n\) = c =const. for any k in the isotropic case. 
The above vectors and matrices are, strictly speaking, of 
infinite dimension, which needs to be truncated for nu- 
merical purposes. In this sense there is no loss of general- 
ity in assuming a smooth /Lt(x) in the course of derivations 
in JFlI Implementation of Eq. (128]) i consists of two steps. 



■ St □ Ep 

c MM num.: 

-JV= 1 
-N = 10 
N = 20 




cpwe num.: 
- * - N = 1 

— o— N = 10 
— N = 20 
(22): 



CMM 



1 / 



FIG. 1: Effective speed c versus concentration / of square rods 
for 2D St/Ep and Ep/St lattices (see details in the text). 



Step 1. Calculate the multiplicative integral (28 )a defin- 
ing Mo. For an arbitrary /i(x), one way is to use a 
discretization scheme. Divide the segment x\ € [0, 1] 



into Ni intervals 



„(0 „(*+!) 



i = l..Ni, of small 
Fourier coefficients 
p„(x^'), n = -N..N and the (2JV+1) x (2JV+1) matrices 



enough length. Calculate 



) = 
2N 



1 



Qo(x^) for each i = l..Ni, and then use the approximate 



formula M, 



o 



m ex P 



A t Q Q {x[ l) ) . Recall that / 



satisfies the chain rule and is exactly equal to exp(AiQo) 
for x\ G Aj if jtx(x) does not depend on x\ within Aj. 
Therefore the calculation is much simpler in the com- 
mon case of a piecewise homogeneous unit cell with only 
a few inclusions of simple shape (see the example below). 
Step 2. Solve the system (M — I)w-j = iwg for unknown 
Wj. First remove one zero row and one zero column in 



the matrix Mq — I (see the remark below ( 10 ) ) . Then the 



vector is uniquely defined and may be found by any 
standard method. Note that only a single component of 
h is needed to evaluate wg_ ■ wj. Finally dividing by (p) 



yields the desired result ( 28 1 



As an example, we calculate the effective shear-wave 
speed c versus the volume fraction / of square rods pe- 
riodically embedded in a matrix material forming a 2D 
square lattice with translations parallel to the inclusion 
edges. A high-contrast pair of materials is chosen such 
as steel (= St, with p = 7.8 10 3 kg/m 3 , p = 80 GPa) and 
epoxy (= Ep, with p = 1.14 10 3 kg/m 3 , p = 1.48 GPa). 
We consider two conjugated St/Ep and Ep/St configu- 
rations, where the matrix and rod materials are either 
St and Ep or Ep and St, respectively. The results are 
displayed in Fig. 1. The curves cmm(/) are computed 



by the present monodromy-matrix (MM) method, Eq. 
( 28 1 1 , they are complemented by the approximation ( 22 1 . 
Also shown for comparison are the curves cpwe(/) com- 
puted from the truncated formula^ of the conventional 
PWE method based on a 2D Fourier transform of (3|. 
Calculations are performed for a different fixed number 
2N + 1 = d of the ID Fourier coefficients of /i(x), which 
implies 2d x 2d monodromy matrix in (28)i and, by con- 
trast, d 2 x d 2 matrix in the PWE formula 4 . Apart from 
this advantage of the MM calculation, it is also seen to 
be remarkably more stable - with a reasonable fit pro- 
vided already at N = 1. The difference between the MM 
and PWE numerical curves is especially notable for the 
case of densely packed steel rods. Interestingly, the MM 
computation and estimate both predict a steep fall for 
c (/) when a small concentration 1 — / of epoxy forms a 
'quasi-insulating network'. The PWE fails to capture this 
important physical feature for reasons described next. 

The far superior stability and accuracy of the MM 
method observed in Fig. 1 can be explained as fol- 
lows. The PWE formulcP implies calculating M\\ f=s 
E| g j< d B s |gf 2 (\g 2 \ + l)' 2 + O {d- 1 ) with bounded co- 
efficients B %1 where g are the 2D reciprocal lattice vec- 
tors (we use here that the components of the vector d±p 
for piecewise constant p(x) are of order (\g2\ + 1) , 
and that the matrix corresponding to Cq 1 is close to 
diagonallv-dominant and hence its eigenvalues are of or- 
der |g| ). Thus the accuracy of the PWE method 
is expected to be of order d^ 1 . In contrast, the accu- 
racy of the MM method, where the ID Fourier expan- 
sion is performed inside a multiplicative integral that is 
'close' to exponential, is expected to be on the order e~ d . 



This can be understood from the MM equation ( 28 1 1 
where the 2d x 2d matrix (M — I) -1 can be replaced by 
2 (Mo — Mq ) _1 with eigenvalues of order e~ n , n = l..d. 
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